%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This Matlab program is for Figure 1 from the Supplemental Appendix of 'The explicit formula for the Hodrick-Prescott filter in finite sample' by Adriana Cornea-Madeira
% HP based on the direct inversion of the matrix (in (2))
clear all; clc;
alpha = 1600; %smoothing parameter
nnn = (100:100:10000)';
Kk = 1; %number of samples
t = zeros(length(nnn),1);
for in=1:length(nnn)
tic
n = nnn(in,1);% sample size

p = n:-1:1;
P = eye(n);
P = P(p,:); % permutation matrix
Q = diag(2*ones(n,1))+diag(-1*ones(n-1,1),1)+diag(-1*ones(n-1,1),-1);% or use tril&tril triu(ones(4,4),-1)
QQ = Q*Q;
g = [-2;1;zeros(n-2,1)];
gg = g*g'; PggP = P*g*g'*P;
E = eye(n);
F1 = E+alpha*QQ;
F = F1-alpha*(gg+PggP); % matrix to be inverted 
Finv = inv(F); % or Finv = F\E;
for k=1:Kk
  % local linear trend model (LLT)
  beta = zeros(n,1);
  tau = zeros(n,1);
  y = zeros(n,1);
  sigma_c = 2;
  sigma_eta = 0; % when sigma_eta = 0, smooth trend model
  sigma_zeta = 0.025*sigma_c; 
  beta(1) = randn(1,1); tau(1) = randn(1,1); y(1) = randn(1,1);
  epsi_1 = randn(n,1); epsi_2 = randn(n,1); epsi_3 = randn(n,1);
  for i = 2:n
      beta(i) = beta(i-1) + sigma_zeta*epsi_1(i,1);
      tau(i) = tau(i-1) + beta(i-1) + sigma_eta*epsi_2(i,1);
      y(i) = tau(i) + sigma_c*epsi_3(i,1);
  end
tau_est = y - Finv*y;  
end
t(in)=toc;
end
figure(1)
plot(nnn,t,'g'); hold on;